Complete high-precision entropic sampling 

Ronald Dickman* and A. G. Cunha-Netto^ 

O ■ Departamento de Fisica, Instituto de Ciencias Exatas, 

(N ■ 

and National Institute of Science and Technology for Complex Systems, 



3 
^1" 



On 



cd 



Universidade Federal de Minas Gerais, C.P. 702, 
30123-970 Belo Horizonte, Minas Gerais, Brazil 



^h ■ Abstract 

O 

Monte Carlo simulations using entropic sampling to estimate the number of configurations of a 
given energy are a valuable alternative to traditional methods. We introduce tomographic entropic 

a 

sampling, a scheme which uses multiple studies, starting from different regions of configuration 
space, to yield precise estimates of the number of configurations over the full range of energies, 
without dividing the latter into subsets or windows. Applied to the Ising model on the square 
lattice, the method yields the critical temperature to an accuracy of about 0.01%, and critical 
exponents to 1% or better. Predictions for systems sizes L = 10 - 160, for the temperature 
of the specific heat maximum, and of the specific heat at the critical temperature, are in very 
close agreement with exact results. For the Ising model on the simple cubic lattice the critical 
temperature is given to within 0.003% of the best available estimate; the exponent ratios fijv and 



r**"- ' 7/1/ are given to within about 0.4% and 1%, respectively, of the literature values. In both two and 

o : 

three dimensions, results for the antiferromagnetic critical point are fully consistent with those of 

the ferromagnetic transition. Application to the lattice gas with nearest-neighbor exclusion on the 

• i— 1 . 

KS ' square lattice again yields the critical chemical potential and exponent ratios [$/v and 7/1/ to good 

5-j ' 



precision. 
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I. INTRODUCTION 

A simple and reliable Monte Carlo method, capable of estimating thermodynamic prop- 
erties at any desired temperature using a single simulation, has been sought for many years 
lrM|. Such a method could in principle offer a great advantage over Metropolis Monte Carlo 
and cluster algorithms, which require a separate study for each temperature of interest. The 
information required to study the full temperature range is contained in the number of 
configurations Cl(E, L) as a function of energy E, and associated microcanonical averages 
[for example, the mean magnetization (m)(E,L)], on a lattice of L d sites. Wang-Landau 
sampling (WLS) [5|, |6| is currently the most widely used such entropic sampling method. 



!T 



Various improvements to the basic WLS scheme have been suggested |7H10|. In WLS of 
larger systems the range of energies (or of other relevant quantities) is divided into subsets 
and each such "window" is sampled separately. The latter procedure has been found to 

rin 

distort the estimates for Q(E,L) in some cases |8l. llOl]. 

Despite the broad interest in WLS, relatively few works have been published in which 
critical exponents are determined using the method. Important exceptions to this general 



trend are the studies of Malakis and coworkers llMl5|. who use a modified WLS procedure 
to find the critical exponents of several pure and disordered spin models. Their method 
involves a two-stage process, first using WLS with the energy range divided into windows 
to identify the critical energy subspace (i.e., the range of energies needed to sample the 



critical region adequately), followed by WLS restricted to the critical range. In Ref. 22 1 
the present authors obtained fair results for the critical exponents of the lattice gas with 
nearest-neighbor exclusion, using WLS with adaptive windows. 

Here we describe a high-precision entropic simulation method that samples the full con- 
figuration space without using windows. Two fundamental aspects of the method are: (1) 
successive refinement of an initial estimate for the Cl(E, L) on the basis of simulations start- 
ing from a diverse set of initial configurations, and (2) use of the final result at size L to 
generate a good preliminary estimate for Cl(E, L') for size V > L. After describing the 
method in Sec. II, we present, in Sec. Ill, applications to the Ising model on the square and 
simple cubic lattices, and to the lattice gas with nearest-neighbor exclusion. A summary 
and conclusions are given in Sec. IV. 



II. METHOD 

Consider a statistical model with a discrete configuration space, and let $ denote a vari- 
able (or set of variables) characterizing each configuration, such as energy or particle number. 
For a given system size, knowledge of the number Cl(i9) of configurations, for all allowed val- 
ues of $, permits one to evaluate the partition function and associated thermal averages for 
any desired temperature. Entropic sampling methods use Monte Carlo simulation to furnish 
estimates of the configuration numbers, which we denote by fi($), reserving f2($) to denote 
the exact values, which are in general unknown. The basic idea of entropic sampling is as 
follows. Let E(C) denote the energy of configuration C. Consider a simulation method that 
generates configurations from the probability distribution 

P(C) oc ■— ^ , (1) 

n[E(c)] 

i.e., the probability of configuration C is inversely proportional to the number of configura- 
tions having the same energy. This could be done using an acceptance probability 



p(C') = min 



n(E) i 



(2) 



n(E> 

(As usual, when the new configuration C is rejected, the current configuration C is counted 
again in the sample.) Let H(E) be the number of times energy E occurs in the sample. In 
a simulation that follows Eq. ([I]), (H(E)) = Const., independent of E. Of course, the Cl(E) 
are in general unknown. If we run the simulation using a set of estimated values, 0,(E), 
in place of Cl(E) in Eq. (OQ), then the better the estimate, the more uniform the histogram 
H(E) should be. The expected number of visits to energy E when we use an estimate Q(E) 
to define p a (C) is 

WE)) cc || (3) 

Suppose we begin using a guess, Qo{E), for Cl(E). Simulating for a sufficiently long time, 
so that H(E) m (H(E)) for all energies E, we may refine our guess, using Eq. ([3]) to write 

Q 1 (E) = 5^W) (4) 



where H i (E) is the number of times energy E occurs in the simulation, and H\ is the average 
of H\(E) over all allowed energies. Equation (j4]) is the basis of our simulation method. The 
idea is to iterate this procedure, starting from an initial guess, enough times so that the 
final approximation, fijv, is very close to the true distribution. (During the final iteration, 
the microcanonical averages required to calculate thermal averages are also determined.) 

For the method to yield a reliable estimate for Cl(E), it is essential that the space of 
configurations be sampled long enough, and broadly enough, that at each step j = 1, ..., N, 
the observed histogram Hj(E) is close to its expected value, that is, Hj(E) oc Q(E)/Qj-i(E), 
for all allowed values of the energy. If this condition is satisfied at each iteration j, we expect 
the estimates Qj^i(E) to approach the exact values. 

The preceding discussion suggests that, to converge, the configuration space must be 
sampled well at each iteration. Our numerical experiments confirm that better results are 
obtained using a relatively small number of iterations, with a relatively large number of 



lattice updates (LUDs) 16|] than with a larger number of iterations, each having a smaller 
number of updates, as shown in the following section. In the studies reported below, for 
example, we use N = 5 steps, each consisting of a set of ten simulations, each of which runs 
for Njj = 10 7 LUDs. 

To ensure adequate sampling of the full range of energies, the simulations in a given 
iteration use a variety of initial configurations. For the Ising model on the square and 
simple cubic lattices, for which the number of configurations is invariant under (1) inversion 
of all spins, and (2) inversion of all spins on one sublattice, we use the following set of initial 
conditions: (1 and 2) random spin orientations; (3) all spins up; (4) all spins down; (5) all 
but two neighboring spins up; (6) all but two neighboring spins down; (7 and 8) all spins 
on one sublattice up, the rest down; (9 and 10) configurations (7) or (8) with a pair of 
neighboring spins inverted. Thus there are four low-energy initial configurations, four of 
high energy (but of low energy for the antiferromagnetic model), and two that have E rs 0. 
Taking a certain liberty with terminology, we call this approach tomographic sampling of the 
distribution Qj, since studies using different starting points are pooled at each iteration. 

A key question concerns the initial estimate, Qq, for the number of configurations. Here 
it is important to note that, in terms of the energy density e = E/L d , the microcanonical 
entropy density, 



s(e,L) = ^- d \nQ(L d e,L) (5) 

depends on system size L rather weakly. Thus the final result for s(e, L) serves as a rather 
good first approximation for a larger system size, V ' . We use a simple initial guess (for 
example, Q(E) uniform over the set of allowed energies, or a mean-field approximation) 
for some small initial system size. For small L the simulation rapidly converges to a good 
estimate of Cl(E,L), and the resulting entropy density s(e,L) is then used to generate the 
initial guess Qq(E,L') for the next system size in the sequence. 

A. Implementation 

All the studies reported here use periodic boundaries; other boundary conditions are 
readily implemented. The Ising spin configurations are classed by the bond number n, i.e., 
the number of nearest-neighbor pairs of spins having the same orientation. On the square 
lattice the allowed values are n = 0,4,6, ...,n max — 4:,n ma x with n max = 2L 2 . (The energy 
is E = 2(L 2 — n).) To begin, we set the initial distribution Q (n). Due to the extremely 
large values these numbers may take (already ~ O(10 30 ) for L = 10), they are represented 
in the form of their natural logarithm. Subsequently, the simulation procedure is called N 
times. (In most of the studies reported here, N = 5.) The simulation uses a single-spin-flip 
dynamics. The acceptance probabilities Qj(n)/Qj(n') (for a transition from energy n to nf), 
are stored in a look-up table. (The allowed values of An are 0, ±2, and ±4.) In addition to 
storing the orientation of each spin, we store the number of neighbors having the same spin 
orientation, as this facilitates evaluation of An. 

Each time the simulation subroutine is called, the histogram Hj(n) is set to zero. The 
latter then accumulates the number of visits to each class n over N sim = 10 simulations (with 
varied initial configurations), each consisting of iVjy LUDs. (In principle, these studies could 
be carried out in parallel, on separate processors.) Once all N sim studies are done we update 
Q(n) as per Eq. (J3J). The acceptance probabilities are recalculated using the new estimate, 
and we proceed to the next iteration, for a total of N iterations. Each time a configuration 
C is generated, be it a new one (following an accepted spin flip) or a repetition (if the new 
configuration is rejected), we update the running sums used to calculate the microcanonical 
averages of the absolute magnetization \M\, as well as M 2 and M 4 . The resulting averages 



are saved, along with \nQ j(n), at the end of each iteration. (Note however that only the 
results from the final iteration, j = N, are used to calculate the thermodynamic properties 
reported below.) 

As mentioned above, the initial distribution f2 (n, L) is generated on the basis of the final 
estimate for a smaller system size. For the smallest size (L = 10 in the case of the square 
lattice) we use a simple approximation for the entropy density, 

S (p)^-ln2(p-l) 2 (6) 

where p = n/L 2 . This is motivated by the fact that Cl(n = L 2 ) ~ 2 L whereas Cl(n = 2L d ) = 
Q(n = 0) = 2, and that s(p) is symmetric about p — 1. (For convenience we subtract In 2 
from the original definition, s = [In Cl(n)] j L 2 .) For L = 10, the distribution after the first 
iteration, f2i(n), is already very close to the final one, showing rapid convergence. 

The quantity h(p) = [H (p) — H}/H characterizes the relative deviation of the histogram 
from a uniform distribution; its evolution at successive iterations is illustrated in Fig. [H 
The result of the first iteration (starting from the parabolic initial approximation) shows 
relatively large variations, while the subsequent histograms are all flat on the scale of this 
figure. The detailed evolution in the subsequent iterations is shown in the inset; the final 
histogram is flat to about 99.5%. 

FIG. 1: (Color online) Ising model, square lattice: relative deviation of histogram, h(p) = [H(p) — 
H]/H, versus bond density p for L = 10. The black curve is the result of the first iteration, while 
the horizontal line corresponds to the four subsequent iterations. Inset: Detail of iterations 2-5. 

For the passage from one system size (L), to the next (L'), we require an algebraic 
approximation for the entropy density s(p,L). Given that Cl(n) = Cl(n max — n), we first 
symmetrize the result about p — 1, and average over pairs of neighboring values to suppress 
small systematic oscillations in s. (The latter are visible in Fig. 2, near p = and p = 2.) 
This defines the symmetrized and (slightly) smoothed entropy density, 

s(p') = \ Hn/L 2 ) + s((n - l)/L 2 ) + s{2 - n/L 2 ) + s(2 - („ - 1)/L 2 )] (7) 

where p' = (2n — 1)/(2L 2 ). For smaller sizes, s is well approximated by a polynomial fit 
in even powers of x = p — 1. For larger systems, however, maintaining a good fit requires 



increasing the order of the polynomial. The fit tends to be worst near x — ±1. (In fact, one 
may anticipate a weak singularity, s ~ |1 — x\ In |1 — x\ at these limits.) We find that, given 
the symmetry s(— x) = s(x), a more convenient alternative is a Fourier cosine series: 



s(x)-s(l) = J2 

j=0 



r (2j + l)irx 



dj cos 



/(*) (8) 



The expansion coefficients are readily calculated using the discrete approximation to the 
standard integral expression; the quality of the fit varies considerably with the number of 
terms in the series. We calculate the maximum absolute deviation 5 = max x \s(x) — f(x)\ 
as a function of the number of terms, and use the value that minimizes 5. For example, for 
L = 120 (square lattice), we obtain the smallest 6 (2.7 x 10~ 5 ) for a series with J = 9151. We 
use the Fourier series for all studies of the two- and three-dimensional Ising models reported 
below. 

The execution time (for a given number of iterations, initial configurations, and lattice 
updates) scales as the number of lattice sites. The largest studies reported here (for system 
size L = 160) required approximately one week on a processor with a speed of 2.8 GHz. 

III. RESULTS 

A. Ising model on the square lattice 

We apply tomographic sampling to the Ising model on the square lattice, for system sizes 
L = 10, 20, 40, 80, 120, and 160. The procedure is as outlined above, except that for the two 
largest sizes, following the five studies using N v = 10 7 lattice updates, a sixth study using 
Njj = 2 x 10 7 is performed. In all cases, the results are calculated on the basis of the final 
study. Means and uncertainties are calculated using five or in some cases six independent 
simulations. 

The number of configurations fl(p) takes its maximum at p = 1 and is (except for 
fluctuations) symmetric about this point. In Fig. |2] we plot the microcanonical entropy 
density 



for several system sizes. At this scale, finite-size effects are only evident near the limiting 
bond densities. 

FIG. 2: (Color online) Microcanonical entropy density s(p) versus p for system sizes L = 10, 20, 
40, 80, and 160 (upper to lower). The inset is a detailed view near p = 2. 

To perform a quantitative test of the method we compare results for the specific heat, mag- 
netization, susceptibilit y, a nd reduced fourth cumulant with known values and with finite- 

]. We define the specific heat per site as c = var (E) / (L d kBT 2 ) , 



size scaling theory 



17 



where E is the total energy. (From here on we employ units such that Boltzmann's constant 
ks = 1.) Figure [3] shows c(T) (for L = 80) as obtained in five independent simulations. 
The curves are indistinguishable on the scale of the main graph, but a detail of the critical 
region does reveal some scatter in the specific heat maxima and the temperatures at which 
they occur. From the analysis of Ferdinand and Fisher [19j one knows that on a square 
lattice of L x L sites (with periodic boundaries), the specific heat takes its maximum value 
at temperature T(c max ) ~ T c (l + 0.3603/L), where T c ~ 2.269185 is the critical temperature. 



Combining the results of Refs. 19| and 20j one also has an expansion for the specific heat 
at the critical temperature: 

c(T c ,L) = A \nL + C + ^ + ^ + 5§ + • ■ ■ (10) 



where A = (2/7r)[ln(l + v / 2)] 2 — 0.494358 and, for an LxL lattice with periodic boundaries, 
C ^ 0.138149, C ± ~ -0.170951, C 2 ~ 0.018861, and C 3 ~ 0.056765. In Fig. |we plot 
c(T c ) versus InL; a least-squares linear fit yields c(T c ) = 0.498(2) InL, i.e., about 1% above 
the exact amplitude. A detailed comparison with theoretical predictions is given in Table 
HI simulation results agree with theory to within uncertainty. Linear extrapolation of the 
simulation results for T(c max ) (for L > 20) versus 1/L yields T c = 2.26966(8), that is, about 



0.02% above the exact value of 2.269185... 2l|. A similar analysis of the energy per site 



at T c yields e c = —1.4147(5), consistent with the exact value, e c = — v2 [2l|. As noted 
above, our sampling method represents the antiferromagnetic (AF) critical point with the 
same precision as the ferromagnetic one. We have verified that T(c max ), c max , and c(T c ) 
associated with the AF transition agree to within uncertainty with the values quoted in 



Table [U for the ferromagnetic transition. 

FIG. 3: (Color online) Specific heat per site c versus temperature for L = 80 in five independent 
simulations. 



FIG. 4: (Color online) Specific heat per site at the critical temperature c c versus InL. The slope 
of the regression line is 0.498(2). 



L 


T(c max ) (th) 


T(c max ) (sim) 


c(T c ) (th) 


c(T c ) (sim) 


10 


2.35 


2.34450(6) 


1.2600 


1.2597(3) 


20 


2.310 


2.30806(8) 


1.6112 


1.6121(5) 


40 


2.290 


2.2889(2) 


1.9582 


1.9573(18) 


80 


2.2794 


2.2793(2) 


2.3031 


2.302(3) 


120 


2.2760 


2.2761(1) 


2.5043 


2.504(5) 


160 


2.2743 


2.2743(2) 


2.64695 


2.655(7) 



TABLE I: Square lattice: comparison of theoretical and simulation results for the temperature of the 
specific heat maximum, and for the specific heat at T c . Note that the theoretical expressions for T(c max ) 
and c(T c ) are subject to corrections of order 1/L 2 and 1/L 4 , respectively. 



We turn now to the magnetization, the susceptibility, and Binder's cumulant. The mag- 
netization per site is given by 



m(T,L) = j^-^^2n(E,L)e-P E (\N + -N4)(E,L) 



(11) 



where /3 = l/Zc^T, Z is the partition function, and (\N + — N_\)(E,L) denotes the micro- 
canonical average of the absolute value of the total magnetization ( JV+ is the number of sites 
with spin a = 1). Results for the magnetization are plotted in Fig. |5j The uncertainties, 
which are not visible on the scale of the main graph, are shown in the inset; as is to be 
expected, they are most severe in the vicinity of the critical point. We note however that 
in all cases, the relative uncertainty in m is less than 0.4%. The finite size scaling (FSS) 
relation m(T c , L) ~ L~^l v permits one to estimate the associated exponent ratio. Using the 
data for all system sizes we find (5/v = 0.1237(5), that is, 1% smaller than the exact value 
of 1/8. 

The susceptibility x = var(M)/(L a! T), and its uncertainty are plotted in Fig. [6J The 
relative uncertainty in x is largest (about 1%) in the critical region, for L = 160; for smaller 



FIG. 5: (Color online) Magnetization per site versus temperature for system sizes L = 10, 20, 40, 
80, and 160. Uncertainties are plotted in the inset. 

systems and away from the critical region, it is considerably smaller. We list simulation re- 
sults for the temperature of the susceptibility maximum, T(j(max), and for the susceptibility 
at T c , in Table HT1 Extrapolating the values of T(x m ax) we find T c = 2.26926(10). Analysis of 
the susceptibility data using the FSS relation x(T c , L) ~ L 7 ^ yields 7/V = 1.754(2), which 
is 0.2% higher than the exact value of 7/4. It is interesting to note that, restricting the anal- 
ysis to system sizes 20 < L < 120, we obtain the slightly superior estimates /3/u = 0.1240(8) 
and 7/1/ = 1.748(2). 



L 


■* \Xmax) 


X(T C ) 


10 


2.4770(4) 


1.7894(4) 


20 


2.3720(1) 


6.093(5) 


40 


2.3208(2) 


20.44(8) 


80 


2.2949(2) 


68.95(25) 


120 


2.2864(1) 


139.5(8) 


160 


2.28196(14) 


233.2(1.8) 



TABLE II: Square lattice: simulation results for the temperature of the susceptibility maximum, and 
for the susceptibility at T c . 



FIG. 6: (Color online) Susceptibility per site versus temperature for system sizes L = 10, 20, 40, 
80, and 160. Uncertainties are plotted in the inset. 



An independent estimate of the critical temperature is afforded by Binder's reduced fourth 
cumulant 23| . Q4 = 1 — (M 4 )/(3(M 2 ) 2 ). Analysis the crossings of Binder's cumulant for 
sucessive pairs of system sizes yields the temperatures and cumulant values listed in Table 
IIHI Linear extrapolation of the results for the last four pairs yields T c = 2.2694(2) and 
Q = 0.610(1). The latter is in good accord with the literature value Q = 0.61071(2) [33 ] . 
Pooling the estimates for the critical temperature derived from the analysis of the specific 
heat, the susceptibility, and the reduced cumulant, we find T c = 2.2695(1), about 0.01% 
above the exact value. 

Our estimate for T c is, as noted, based on extrapolations of three sets of size-dependent 
pseudocritical temperatures, i.e., T(c max ), T(xmax), and T x . These quantities are plotted 



versus 1/L in Fig. [7J We analyzed the behavior of the pseudocritical temperature associated 
with the susceptibility maximum, which should follow T(xmax) — T c ~ L' 1 ^. Using our 
estimate of T c = 2.2695, a linear fit of ln[T(x max ) — T c ] versus InL yields v = 0.99(1) 
consistent with the critical exponent v — 1 for the two-dimensional Ising model. 

FIG. 7: (Color online) Two-dimensional Ising model: pseudocritical temperatures (upper to lower) 
T(xmax), T(c max ), and T x versus inverse system size. Error bars are smaller than the symbols. 



L, V 


T x 


Qx 


10, 20 

20,40 

40, 80 

80, 120 

120, 160 


2.2635(1) 
2.2692(2) 
2.2689(4) 
2.2696(4) 
2.2693(5) 


0.61372(10) 

0.61115(20) 

0.6115(10) 

0.6096(17) 

0.611(2) 



TABLE III: Ising model, square lattice: simulation results for the temperature and the cumulant value 
at crossings of Binder's cumulant for pairs of successive system sizes. 

As noted above, we find it preferable to perform a relatively small number of long itera- 
tions than a large number of shorter ones. To illustrate this point, we performed a series of 
studies using a total of 5 x 10 8 LUDS, divided into 5, 10, 25 and 50 iterations, for system size 
L = 20. In Fig. [8] we plot the relative uncertainties (estimated using the standard deviation 
calculated over a set of five independent studies), of T(c max ), T(xmax), and of the magne- 
tization, susceptibility, and Binder's cumulant, evaluated at the critical temperature. As 
is clear from the figure, the uncertainties decrease systematically as we reduce the number 
of iterations; the relative uncertainty is roughly 5-10 times larger using 50 iterations as 
compared to only five. The inset, for m c , illustrates the general trend: despite the mono- 
tonic variation of the uncertainties, the estimates obtained using different values of N are 
mutually consistent. 



FIG. 8: (Color online) Relative uncertainties in T(c max ) (diamonds), T(xmax) (x), m c (filled 
squares), Xc (open squares), and Q c (circles) for L = 20, versus number of iterations N, in studies 
using a total of N x Njj = 5 x 10 7 lattice updates. Inset: estimates for m c versus N . 



B. Ising model on the simple cubic lattice 

The procedure is essentially the same as for the square lattice. We begin sampling at 

L — 4, using the simple initial guess for s(p) given in Eq. ([6]). Subsequently we study system 

sizes L = 8, 12, 16, 20, 24, 28, 32, and 36. The smaller increment in L (as compared with 

the two-dimensional case), is motivated by faster convergence, and by the need to study a 

reasonably large number of system sizes in the FSS analysis, which in this case includes a 

correction to scaling term. 

FIG. 9: (Color online) Ising model, simple cubic lattice: specific heat per site versus temperature 

for system sizes L = 8, 12, 16, 20, 24, 28, 32, and 36. 

FIG. 10: (Color online) Ising model, simple cubic lattice: pseudocritical temperatures T{c max ) 
(lower) and T(x m ax) (upper) versus \jL x ' v . Inset: temperatures T x of Binder cumulant crossings 
versus 1/L (here L denotes the geometric mean of the system sizes). 

The specific heat in the critical region is shown in Fig. [9j the relative uncertainty in 
these results is at most 0.5%. We estimate the critical temperature on the basis of T(c max ) 
and T(xmax) (see Fig. [10]). Using the values of T(c max ) for L = 12 - 36, a quadratic 
least-squares fit versus 1/L 1 ^ yields T c = 4.5116(1); a similar analysis using T(xmax) yields 



4.5114(1). (We use the literature value, v = 0.6301(8) 24J.) These results are in good accord 



with the best estimate of T c = 4.511528(6) 25] . The energy per site at the critical point, 



e c (L), is essentially linear when plotted versus 1/L 1 /"; extrapolation to infinite size yields 



e c = —0.9928(2), which compares well with the series-expansion estimate, e c = —0.9922 26]. 



For each system size, the values of c max and T(c max ) associated with the antiferromagnetic 

critical point are fully consistent with the corresponding ferromagnetic values. 

FIG. 11: (Color online) Ising model, simple cubic lattice: magnetization (lower), specific heat 

(middle), and susceptibility (upper) evaluated at the critical temperature versus system size. Error 

bars are smaller than symbols. 

The critical specific heat, susceptibility, and magnetization are plotted versus system size 
(on log scales) in Fig. [TTJ In light of the slight curvature evident in these plots, and given the 
smaller system sizes used in the three-dimensional studies, we found it useful to include a 
correction to scaling in the FSS analysis. For example, we fit the data for the magnetization 
at T c using 

In m(T c , L) = --In L + cL~ Vl + Const. = f L + Const. (12) 



where y^ the dominant irrelevant scaling exponent, is taken as 0.8. (Note that the constant 
term is not a fitting parameter, as we determine the values of (3/v and c by minimizing the 
variance between the simulation data and //,.) The best-fit parameters for system sizes L = 8 
- 36 are c = —0.064 and (3/v = 0.521(12). The uncertainty a-s in j3 = /3/v is calculated as 
follows. Let Yi = In m(T c , L) — f^, and let a\ denote the variance of the Yl, for parameters 
/3 and c chosen to minimize this variance. Further, let a\ denote the statistical uncertainty 
in the simulation result for In m(T c , L), and define a\ = a 2 L + a\. Then we take 



T- = V^ 7f 2 ( — 

~l \d^m{T c ,L) 



•4=Elsrans ■ < 13 > 



where the derivatives are determined numerically. The uncertainties in the other exponent 
ratios are estimated in a similar manner. 

For the susceptibility, we expect the dominant correction to FSS to be oc L~ m with 



y 2 ^ 1.96 [24j, and write 



In x(T c , L) = -\nL + d L~ m + Const. (14) 

In this case we find j/v = 1.987(4). (Using a correction to scaling term oc 1/L we instead 
obtain 7/1/ = 1.955(11).) A similar analysis of the specific heat data (for system sizes L = 12 
- 36, using a correction term oc L~ a / U ) yields a/v = 0.161(3). Compared with the literature 



values |25|, fi/v = 0.519(2), 7/z/ = 1.963(3), and a/v = 0.174(4), our results for exhibit 
errors of 0.4%, 1.2%, and 7.5%, respectively. (We note, however, that determination of a/v 
via simulation is in general a difficult task.) 

Extrapolation of the crossings of Binder's cumulant (see Fig. [101 inset), yields T c = 
4.5124(16), which, while less precise than the estimates cited above, is consistent with the 
best estimate for the critical temperature. The asymptotic value of Binder's cumulant 
obtained from the crossings is Q4 = 0.47(1), while extrapolation of the cumulant values 
Q^(T C ,L) yields 0.467(1). The reference value for the three-dimensional Ising model is 



Q 4 = 0.465(3) [24 1 . 



Entropic sampling is particularly advantageous for data collapse analyses, as it furnishes 
thermodynamic quantities as continuous functions of temperature. FSS predicts that mag- 
netization curves for diverse system sizes should fall on universal curves (for T > T c and 
T < T c ) when plotted in the form L /3 / iy m(L 1 / i/ |T — T c \). Figure H~2l shows a near-perfect col- 



lapse of the magnetization data for eight system sizes (L = 8-36); the associated exponents 

are v = 0.62 and /3/v = 0.512. Figure [TBI is a similar plot for the susceptibility, again using 

v = 0.62, and j/u = 1.99. In this case the data for T > T c collapse perfectly while those for 

T < T c approach a common scaling function with increasing system size. (The downward 

curvature in these plots signals the low-temperature boundary of the scaling region.) 
FIG. 12: (Color online) Ising model, simple cubic lattice: scaled magnetization m* = L°' u m 
versus scaled temperature t* = L ' U \T — T c \ for eight system sizes (L = 8-36), using v = 0.62 and 
p /u = 0.512. 

FIG. 13: (Color online) Ising model, simple cubic lattice: scaled susceptibility x* = -^~ X versus 
scaled temperature t* = L l ' u \T — T c \ for eight system sizes (L = 8-36), using v = 0.62 and 
j/v = 1.99. 



C. Lattice gas with nearest-neighbor exclusion 



We apply our entropic 



on the square lattice 



22 



samp 



28 



ing to the lattice gas with nearest neighbor exclusion (NNE) 
29] , using entropic sampling to estimate the number of distinct 
configurations Cl(J\f, L) with M particles satisfying the NNE condition, on an L x L lattice 
with periodic boundaries. The grand partition function is 



Jvrnax 

E(z,L)= ^Tz^A^L), 

A/"=0 



(15) 



where z = e M is the fugacity, /x = J2/k B T (p denotes the chemical potential), and J\f max is 
the maximum possible number of particles, equal to L 2 /2 on the square lattice. (In what 
follows we refer to \x as the chemical potential.) The model exhibits an Ising-like phase 
transition on bipartite lattices; in the ordered phase a majority of the particles occupy one 
of the sublattices. The associated order parameter is the difference between the occupancies 
of sublattices A and B: 



M ir 



o\ 



(16) 



xeA xeB 

where <r x is the indicator variable for occupation of site x. 

We begin our study with a small system (L = 8), using a flat initial distribution, In Q = 0. 
Following the procedure described in Sec. II, we perform N = 5 iterations, with increasing 
numbers of lattice updates [30]. We use ten initial configurations at each iteration: five 



corresponding to an empty lattice, and five in which one sublattice is fully occupied. Here 
each trial move is either an insertion or a removal of a particle; target sites are chosen at 
random. For system sizes L' > 8 we construct the initial guess for Cl(J\f, L') using a tenth- 
degree polynomial fit to s(p,L) = [\nCl(Af, L)]/L 2 , where L is the previous size studied and 
p = N I ' L 2 . We study six sizes in the range 8 < L < 120. 

Figure dH shows the susceptibility, x(a0 = -^ 2 ((0 2 )« ~~ (0)u)> as a function of the chem- 
ical potential for different system sizes; the inset is a similar plot of the order parameter. 
Estimates for the critical chemical potential p c are obtained via analysis of the maxima of 
the susceptibility and of the compressibility, n(p) = L 2 ({p 2 ) — (/?)„)/ (/>)„• Extrapolation 
versus 1/L yields p CtX = 1.3359(3) and /i C;K = 1.344(2). The former, more precise value, 
is about 0.14% higher than the literature value. FSS analysis of the susceptibility maxima 
yields 7/V = 1.750(2), consistent with the exact value for the two-dimensional Ising model. 

Using the high-precision result for the critical chemical potential obtained by Guo and 



B16te [29j, p c = 1.33401510027774(1), we calculate pJL), </> c (L) and Q C (L), where Q 



J ) / (0 4 ) 1S related to Binder's reduced cumulant 23[. Linear extrapolation of these 

and_9 c = 0.8565(5). These 
. (Using our own less 



29 



quantities versus 1/L yields p c = 0.36773(1), f3/v = 0.1247(3' r 

results are in very good agreement with the literature values 

accurate estimate, p c = 1.3359(3) we obtain f3/v = 0.122(2) and p c = 0.36815(3).) We 

summarize our main results for the lattice gas and compare them with reference values in 

Table IV. 

FIG. 14: (Color online) NNE lattice gas: susceptibility versus chemical potential. Inset: order 
parameter versus chemical potential, system sizes as indicated. Error bars are smaller than the 
symbols. 



D. Critical singularity of the entropy density 

The critical point corresponds to a value of p such that, in the thermodynamic limit, 
the second derivative s" = d 2 s/dp 2 = 0. Since direct estimation of the second derivative of 
numerical data via finite differences is not a viable procedure, we instead use a Gaussian 
filter, 





ES 


AWWLS [22; 


Literature values 


^c,\ 


1.3359(3) 


1.330(1) 


1.33401510027774(1)° 


r^c,K 


1.344(2) 


1.337(2) 




Qc 


0.8565(5) e 


0.852(6) 


0.856°; 0.855(l) c ; 0.85625(5)°' 


Pc 


0.36773(l) e 


0.36800(5) e 


0.3677429990410(3)° 


l/v 


1.750(2) 


1.762(8) 


7/4 (exact) 


0/u 


0.1247(3) e 


0.123(2) e 


1/8 (exact) 



°Guo and Bldtc|2 
6 Burkhardt and Dcrri 
c Nicolaides and Bruce 
d Kamieniarz and B16tel3 




6 Values obtained using \x c from [29 1. 
Ajsing our best value for fi c we obtain fl/v 
P/v = 0.130(9) and p c = 0.3675(5). 



= 0.122(2), p c = 0.36815(3) while the AW study yields 



TABLE IV: Critical values for the NNE lattice gas obtained via Entropic Sampling (ES) and WLS 
with adaptive windows (AW). The results of 32| were obtained via Monte Carlo simulation, those of 



3jJ, |33j) via transfer-matrix analysis. 



s"(p;a) 



dp's(p')g"(p - p'; a) 



(17) 



where g(p; a) is a normalized Gaussian distribution with mean zero and standard deviation 
o. (The filter width a is chosen large enough to suppress fluctuations but small enough 
to maintain resolution; it ranges from about 5Ap for L = 20 to about 50Ap for L = 160, 
where Ap = 1/(2L 2 ) is the bond density increment.) Figure [T5l for the square lattice, shows 
that \s"\ exhibits a sharp minimum, which appears to approach zero as L increases. The 
minimum near p — 1.7 corresponds to the ferromagnetic transition; note that the expected 
(infinite-size) value is rho c = 1 — e c /2 = 1 + l/v2 = 1.707 107.... The minimum near p = 0.3 
corresponds to the antiferromagnetic critical point. Figures [16] and H7] are a similar plots 
for the Ising model on the simple cubic lattice (for which p c = (3 — e c )/2 ~ 1.997), and the 
NNE lattice gas, respectively. 

FIG. 15: (Color online) Ising model, square lattice: s"{p) versus p for system sizes L = 20, 40, 80, 
120, and 160. The minimum in \s"\ grows sharper with increasing system size. 



FIG. 16: (Color online) Ising model, simple cubic lattice: s"(p) versus p = n/L 3 for system sizes 
L = 16, 24, and 32. The three curves are indistinguishable at this scale. 



FIG. 17: (Color online) NNE lattice gas, square lattice: s"{p) versus p = N/L 2 for system sizes 
L = 20, 40, 80, and 120. 

E. A note on restricted sampling 

Our method, as noted, samples the entire range of energies, and while this is advantageous 
in some circumstances, one might inquire whether the sampling might be restricted or at 
least concentrated in a region of particular interest, for example, to those energies that occur 
with a significant probability in the critical region. In Wang-Landau sampling, imposing 
fixed limits on the sampling range has been found to distort the estimates for Q(E) in 



certain cases 



10 



22| . We therefore consider a smooth preference in sampling, in which the 



acceptance probability, Eq. ([2]) is modified so: 



p a (C) = min 



n(E')f(E>y J 



(18) 



Now the probability of visiting energy level E is oc 1/ f(E), and the formula for updating 
Q(E) becomes, 

n j+1 (E) = ?fflf(E)n j (E). (19) 

In principle, one can penalize visits to "uninteresting" regions using larger values of /, which 
are then compensated for in the modified updating formula. Since a greater fraction of the 
sampling occurs in the region of interest, one expects to obtain better statistics there, for 
the same total number of lattice updates. 

We set / = 1 in the region of interest, that is, for n > n m i n . In one series of studies, 
/ = exp[a(n m in — n)} for n < n m i n , while in the second we set / = exp[a(n m j n — n) 2 ]. Here a 
is a parameter that controls the sharpness of the cutoff in sampling. In the linear case, using 
n min ~ L 2 (i.e., half the maximum bond number), and a in the range 0.1 - 0.5, we obtain 
distributions Q(n) (for L = 40 and 80), that differ significantly from those obtained with 
uniform sampling, even on the region with / = 1. In the second set of studies we chose a 
such that f{n = 0) = e~ 5 ~ 0.007; / grows quite slowly in the vicinity of n min . In this case 
the estimated number of configurations does not exhibit obvious distortions, but we find 
that the susceptibility and specific heat, evaluated near T c , are significantly greater than the 



values obtained with uniform sampling. The relative differences grow with the system size. 
For example, the maximum susceptibility Xmax is found to be 0.6% (1%) higher than that 
obtained using unrestricted sampling, for L = 40 (80). This occurs despite the fact that at 
T c , the the probability e~^ E Q(E) associated with bond numbers n < n min is < 10~ 30 , so that 
the contribution to thermal averages due to this range of n values is completely negligible. 
We conclude that restricting sampling in this manner does not improve the quality of the 
results, and in fact appears to cause systematic errors in thermal averages. 

IV. SUMMARY 

We have devised an entropic sampling method that covers the entire range of energies, 
and yields results of good precision for a modest expenditure of computer time. For the 
two-dimensional Ising model, we obtain the critical temperature to within about 0.01%, 
the exponent ratios 7/V and /3/u to within 1%, and excellent agreement with theoretical 
results for the position and value of the specific heat and susceptibility maxima. Very good 
results for the Ising model on the simple cubic lattice, and for the lattice gas with nearest- 
neighbor exclusion, are also obtained. The method is relatively simple, and avoids sampling 
in "windows" with the attendant problems of patching together the results, and of possible 
distortions at the boundaries. 

The sizes used in the present study are relatively modest, and indeed very large sizes 
are not required to obtain precise results for the Ising model. (One benchmark study of 



the three-dimensional Ising model 24], for example, uses L < 40. As the latter work 
illustrates, precise results for a large set of modest system sizes, coupled with finite-size 
scaling analysis, may be just as useful as simulations of much larger systems.) In the 
present study we observe a reduction in quality (as reflected in the relative uncertainty of 
thermodynamic variables), with increasing system size. We believe that an enriched set of 
initial configurations, permitting more thorough sampling of configuration space, will help 
to maintain precision as one increases the system size. This is an important issue for future 
work. 

Key features of our tomographic sampling method are: pooling results starting from very 
different initial configurations before updating the estimate for Cl, and using the results 
for a given system size to generate an initial estimate for the next larger size. The method 



furnishes globally accurate results, as well as precise estimates for critical properties, using a 
single set of parameter-free simulations, thereby opening new possibilities for the application 
of entropic sampling. 
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